Ionizing terahertz waves with 260 MV/cm from scalable optical rectification

Terahertz (THz) waves, known as non-ionizing radiation owing to their low photon energies, can actually ionize atoms and molecules when a sufficiently large number of THz photons are concentrated in time and space. Here, we demonstrate the generation of ionizing, multicycle, 15-THz waves emitted from large-area lithium niobate crystals via phase-matched optical rectification of 150-terawatt laser pulses. A complete characterization of the generated THz waves in energy, pulse duration, and focal spot size shows that the field strength can reach up to 260 megavolts per centimeter. In particular, a single-shot THz interferometer is employed to measure the THz pulse duration and spectrum with complementary numerical simulations. Such intense THz pulses are irradiated onto various solid targets to demonstrate THz-induced tunneling ionization and plasma formation. This study also discusses the potential of nonperturbative THz-driven ionization in gases, which will open up new opportunities, including nonlinear and relativistic THz physics in plasma.

OR in lithium niobate (LN) has also been widely used due to its high nonlinearity 30 and high damage threshold 31 , enabling the production of high-energy THz radiation.In addition, a tilted pulse front scheme can be adopted to achieve phase (or velocity) matching in a prism-shaped LN 14,30,32,33 .In this scheme, the generated THz wave propagates at the same velocity as the tilted laser pulse front inside the LN so that THz energy can continuously grow with the propagation distance.Recently, THz energy as high as 13.9 mJ has been produced by using two stacked cryogenically cooled LN prisms 33 .However, the radiation frequency produced from an LN prism is mostly peaked below 0.5 THz, which naturally leads to large focal spot sizes ( > mm), consequently limiting the peak THz field strength to at most 7.5 MV cm −1 until now 33 .
Recently, a new phase-matching condition has been discovered in LN, which does not require any pulse front tilting 34 .The working principle is that the velocity of THz waves is generally frequency-dependent and varies so large between two phonon resonance frequencies (7.4 THz and 18.8 THz in LN) that there exists a frequency at which the THz wave and the laser pulse driving OR propagate at the same velocity 34,35 .For an 800-nm optical pulse, this occurs at approximately 15 THz, where the refractive index of LN n THz becomes equal to the optical group index of refraction n g (n THz = n g = 2.3).The resulting phase-matched radiation at 15 THz, in principle, can be tightly focused due to its relatively short wavelength of 20 μm, potentially producing strong electromagnetic fields at the focus.This scheme uses a planar LN crystal and is favorable for scaling up the output THz energy 36 .
In this study, we demonstrate high-energy, highintensity THz generation from large-area planar LN crystals by utilizing phase-matched OR at 15 THz.We also characterize the generated THz waves by their energy, beam profile, spectrum, and pulse duration.Such measurements allow us to determine the peak intensity or field strength of the THz waves at the focus.We also perform numerical simulations on THz wave generation in LN to accurately reproduce the THz pulse shape and spectrum, as well as to understand what thickness of LN can yield the highest THz energy.Finally, we experimentally demonstrate THz-induced ionization by focusing high-power 15-THz waves onto various solid materials, including metals, semiconductors, and polymers.

Experiment
A schematic of our setup is shown in Fig. 2. In this experiment, we used a 150-TW Ti:sapphire laser capable of producing <3.4J, >25 fs pulses at a repetition rate of 5 Hz 37 .The laser was linearly polarized (horizontal) and had a center wavelength of 800 nm with a 78-nm full width at half-maximum (FWHM).The laser pulses after the pulse compressor were directed to the experimental vacuum chamber (see Fig. 2) through a series of vacuum tubes to avoid air-induced self-nonlinearities.The laser was operated in a single-shot mode to minimize cumulative degradation of the compressor gratings due to thermal loading.
Inside the vacuum chamber, the laser beam was apertured by a 60-mm-diameter iris diaphragm and then incident onto a 3-inch-diameter LN crystal wafer for THz generation via OR.The extraordinary axis of the LN wafer was aligned parallel to the laser polarization for maximal THz generation.The LN wafer was slightly tilted to prevent laser back reflection toward the laser system.After the LN wafer, the THz waves were

OR in planar LN
OR in DSTMS [16]   OR in LN Prism [33]   Plasma in solid [24]   CTR [25]   PCA [28]   STE [29]   DFG [12]   Plasma in liquid [23]   Plasma in gas [21]   OR in ZnTe [15]   OR in DAST [17]   OR in OH1 [16]   Plasma in liquid decoupled from the IR laser by a pair of windows coated with indium tin oxide (ITO) layers (250-nm thickness) on the front reflection surfaces.This ITO decoupler was tilted at a Brewster angle of 56°, providing minimal reflectance of <0.0001% at 800 nm and maximal reflection of 64% at 15 THz 34 .For characterization, the generated THz waves were intercepted by a metallic 45°m irror and directed out of the vacuum chamber through a potassium bromide (KBr) window of 3-inch diameter and 25-mm thickness and then focused by a goldcoated, off-axis parabolic mirror (OAP) of 3-inch diameter and 4-inch focal length.More details on THz characterization are provided in the Materials and methods.

THz energy
For THz generation, large-area x-cut congruent LN wafers were used.The LN wafers had a diameter of 3 inches and thicknesses of 25 μm, 50 μm, and 100 μm, either doped or undoped with 5% MgO.The THz energies generated from the LN wafers were measured using a pyroelectric detector, as illustrated in Fig. 3.
Figure 3a, b plot the output THz energy generated from the MgO-doped and undoped LN wafers, respectively, as a function of the group delay dispersion (GDD) of the incident laser pulse.Here, the incident laser energy on the LN wafers was fixed at 2.3 J except the scan marked with green scatters (3.4 J).The voltages measured by the pyroelectric detector were converted to THz energies using the detector's responsivity of 0.255 μJ V −1 at 15 THz.This responsivity was obtained by using a femtosecond 800-nm laser as a calibration source, together with relative frequency-dependent correction data provided by the detector vendor (see Figs. S1 and S2 in the Supplementary Material).The THz energy is retrieved after taking into account the transmission of the filters and attenuators used in the beam path, including the ITO decoupler (64% transmission), KBr window (83%), 7.3-μm longpass filter (49%), 15.7-THz bandpass (42%), 10-mmthick Si window (14.4%) and Si wafers (50.6, 58.8, and 54.1%).Note that the THz energy plotted in Fig. 3 represents the energy estimated immediately after the LN emitter.
Figure 3a, b also show that the output THz energy peaks near zero GDD, indicating that the shorter laser pulse produces more THz energy at a fixed laser energy.Additionally, with increasing LN thickness, the optimal GDD shifts to a more negative value.This is due to a combination of the intrinsic positive dispersion (430 fs 2 mm −1 at 800 nm) in bulk LN and nonzero third-order dispersion (TOD) present in the laser pulses.In addition, Fig. 3a, b show that a thinner LN wafer produces more THz energy.This is mostly due to THz absorption in LN (see the simulation section for details).In addition, MgO-doped LN wafers yield approximately 1.4 ~1.6 times more THz energy than undoped ones for the same thickness.This is because MgO doping can enhance not only the damage threshold 31 but also the nonlinear coefficients d 33 and d 31 38,39 .The green scatter at zero GDD represents the highest THz energy, reaching up to 0.91 mJ, obtained with a 25-μm-thick MgO-doped LN wafer at a pump energy of 3.4 J.
Figure 3c, d show how the output THz scales with increasing laser energy for three thicknesses of MgOdoped and undoped LN wafers.In our experiment, the laser energy was varied by rotating a half-wave plate placed before the laser pulse compressor, with the optimal GDD value fixed for each LN.Since OR is a second-order nonlinear effect based on χ (2) , the output THz energy is expected to scale quadratically with the input laser energy.This quadratic dependence is observed at relatively low laser power, but it slowly deviates with increasing laser power, as shown in the log-log (inset) plots.The signal even saturates with the 100-μm-thick LN wafer.This deviation and saturation is due to three-laser-photon absorption and the cascading effect in LN 34,35,40,41 .The cascading effect broadens the pump laser spectrum and can weaken the OR process, as detailed in the simulation section.

THz beam profile
For spatial characterization, the generated THz waves were focused by the OAP mirror onto a microbolometer camera (see Fig. 2).A 7.3-μm longpass filter and a 15-THz bandpass filter were inserted in the beam path to block any unwanted waves.One or more Si wafers were also used to attenuate the THz waves and prevent possible THz-induced damage to the camera.The LN wafer used here was a 25-μm-thick MgO-doped wafer.Figure 4a presents the evolution of the THz beam profiles captured along the propagation direction at an incremental distance of 1 mm from the focus (z = 0 mm).The beam diameter decreases almost linearly as the camera's position approaches the focal point while maintaining its circular shape.Locally bright edges, however, appear in nearly all profiles, which are attributed to a nonuniform intensity distribution of our laser pulses (see Fig. 4a and Fig. S3).
Figure 4b shows the THz beam profile at the focus.The horizontal and vertical lineouts are well fitted with Gaussian profiles, with focal sizes of w h = 42 ± 1 μm and w v = 44 ± 1 μm in FWHM (see Fig. S4 for fitting results).The geometric average is w FWHM = ffiffiffiffiffiffiffiffiffiffi ffi w h w v p = 43 ± 1 μm in FWHM or w = 37 ± 1 μm (radius at 1/e 2 ).This is ~1.7 times larger than the diffraction-limited value of w = 22 μm, estimated from f/# = 1.7 and the central wavelength of 20 μm.

THz temporal and spectral profiles
To obtain the temporal profile of the generated THz radiation, we utilized a single-shot THz interferometer 42 .In our single-shot scheme, the interference between E 1 ðtÞ and E 2 ðt À τÞ, representing two THz fields separated by τ in time, yields a time-integrated intensity given by 42 IðτÞ ¼ where ð Þ represents a cross-correlation function that depends on τ.The Fourier transform of the crosscorrelation function can be expressed in the frequency domain as which provides the spectral information ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi I 1 ðωÞI 2 ðωÞ p of the THz waves.For two identical THz waves, EðtÞ and Eðt À τÞ, separated by τ in time, the Fourier transform of the autocorrelation term provides the spectral power IðωÞ.
In our geometry, two identical THz beams are incident on the microbolometer camera with a crossing angle of θ (see Fig. 2).The relative time delay τ can be mapped onto the microbolometer's sensor surface as where x is the distance parallel to the sensor surface and c is the speed of light.The crossing angle θ can be deduced from the mapping ratio τ=x, which can be experimentally measured by scanning the delay stage.In our setup, the entire interference fringes were translated by 35 pixels when a 400-fs time delay was imposed.This gives τ=x= 11:4 fs pixel −1 and θ=11.6°.THz intensity profiles, I 1 x; y ð Þ and I 2 x; y ð Þ.Both images show no interference patterns, as expected.Figure 5d represents a noise intensity I 0 x; y ð Þ captured with all beams blocked, representing the background noise.The autocorrelation of the THz field can be obtained from ΔI x; y ð Þ ¼ I x; y ð ÞÀI bg x; y ð Þ, where I bg x; y ð Þ is the background signal corresponding to I 1 þ I 2 in Eq. (1).To take into account shot-to-shot variations in the camera's counts, we model the background signal as , in the time domain.Clearly, the autocorrelation exhibits multicycle oscillations.The corresponding power spectrum is obtained by a Fourier transformation and plotted in Fig. 5g (red line).The spectrum peaks at approximately 14.9 THz with a 1.9-THz FWHM.This bandwidth gives a Gaussian transform-limited pulse duration of 0.23 ps in FWHM.
Alternatively, the spectrum can be directly obtained from I x ð Þ without involving I bg x ð Þ.The black line in Fig. 5g represents the Fourier transform of I τ ð Þ (black line) in Fig. 5f.This yields almost the same spectrum, except for the large low-frequency components due to an uncompensated background signal.Such a low-frequency region (<3 THz) is well distinguished from the frequency band of interest and can be ignored.

Simulation of THz generation
To simulate THz generation via OR, we solved onedimensional (1-D) coupled forward Maxwell equations (FMEs), as detailed in the Materials and methods.First, the refractive index n THz , the linear absorption coefficient α T , and the effective nonlinear coefficient d eff of LN are calculated as a function of THz frequencies 35 and plotted in Fig. 6a.The dotted horizontal line indicates the group index of refraction n g = 2.3 of the laser at 800 nm.Phase matching (n THz = n g = 2.3) occurs at three frequencies, 8.3 THz, 14.8 THz, and 19.3 THz, but only 14.8 THz is expected to be dominant because of its lowest absorption coefficient.d eff is also relatively large at 14.8 THz (d eff ≈ 81 pm V −1 ) to produce high-energy THz radiation.
Figure 6b shows the calculation of the output THz spectral power with an increasing propagation distance within LN.The simulation was performed by setting the input laser pulse to have an energy of 2.3 J, a fluence of 0.08 J cm −2 , a Gaussian spectrum with a FWHM of 78 nm centered at 800 nm, and initial GDD and TOD values of 0 fs 2 and 1000 fs 3 , respectively.It shows that a noticeable  This trend is more evident in Fig. 6c, which plots the cumulative THz energy at 10-20 THz as a function of the initial laser GDD and the propagation distance in LN.The result shows that 15-THz radiation can be most efficiently produced from a thin (<25 μm) LN with the input laser pulse duration as short as possible, i.e., GDD ≈ 0 fs 2 .Figure 6c also plots the expected THz conversion efficiencies (solid lines).As the propagation distance increases in LN, the output THz energy initially rises, reaching its peak at 18 ~25 μm, and then either decreases or saturates depending on whether the cascading and self-phase modulation (SPM) effects (see the simulation section in Materials and methods) are included (red line) or not (yellow line).Together with dispersion, the cascading effect can broaden the laser pulse duration, thereby weakening the OR process.The SPM effect, however, is negligible in our simulation due to its small nonlinear refractive index in LN 34 .Even without these effects, the THz yield saturates because of linear absorption (~0.14 μm −1 at 15 THz), which fundamentally limits the conversion efficiency.The scatters in Fig. 6c are our experimental data obtained with three thicknesses of LN (25, 50, and 100 μm), with (cyan circles) and without (amber squares) MgO doping.The data agree reasonably well with the simulation results.
Figure 6d shows the simulated THz waveform after a 25-μm propagation in LN.It exhibits multicycle 15-THz oscillations on top of low-frequency oscillations at <5 THz 35 .The waveform, extracted after 12 ~18 THz bandpass filtering, is shown in the top inset.The corresponding intensity envelope (red line) is shown in the bottom inset and provides a pulse duration of ~0.23 ps in FWHM. Figure 6e presents the autocorrelation trace (red line) of the simulated THz waveform, together with the measured waveform (black line) for comparison.Their corresponding power spectra are obtained by Fourier transformations and shown in Fig. 6f.The simulated spectrum peaks at 14.7 THz with a bandwidth of 1.4 THz in FWHM, reasonably consistent with the measured peak at 14.9 THz with an FWHM of 1.9 THz.The bandwidths show a 36% discrepancy, which is possibly due to the uncertainty of many parameter values used in the simulation.

THz peak intensity and field strength
With an assumption of a Gaussian spatial and temporal pulse shape, the peak THz intensity is given by 43 , where E is the THz pulse energy, w is the THz beam waist at 1/e 2 , and τ T is the THz pulse duration in FWHM.The highest THz energy achieved immediately after the LN wafer is 0.91 mJ.The energy drops to 0.47 mJ after the ITO decoupler (64% transmission) and the Brewsterangled Si wafer (81%) when focused inside the vacuum chamber.With E = 0.47 ± 0.05 mJ, w = 37 ± 1 μm, and τ = 0.23 ± 0.02 ps at the focus, the peak intensity of the focused THz waves is estimated to be I T = (8.9± 1.3) × 10 13 W cm −2 .The corresponding peak electric and magnetic field strengths are E T = 260 ± 20 MV cm −1 and B T = 87 ± 7 T, respectively (see Fig. S4 for detailed calculations of the uncertainties).To the best of our knowledge, this is the highest peak intensity or field strength demonstrated so far at 0.1 ~20 THz, as shown in Fig. 1.
An intense (>10 13 W cm −2 ) THz pulse, when focused into a neutral gaseous or solid medium, can ionize the constituent atoms or molecules and convert the medium into a plasma.In the plasma, the liberated electron acquires an average kinetic energy in the oscillating THz fields, i.e., the ponderomotive energy , where v o is the quiver velocity of the electron, and m e and e are the electron mass and charge, respectively.In practical units, U p [eV] = 9.33 × 10 −14 I T [Wcm −2 ]λ 2 [μm].At I T = 8.9 × 10 13 W cm −2 , the cycleaveraged kinetic energy of the electron can be as large as U p = 3.2 keV.This energy is more than 5 orders of magnitude higher than the THz photon energy of 62 meV at 20 μm owing to the Iλ 2 -dependent scaling law.Such energy can be converted to produce multi-keV photons by high harmonic generation when the liberated electron recombines with its parent ion 44 .In addition, the normalized THz vector potential, a 0 ¼ eE T = m e cΩ ð Þ¼v=c, gets close to 1 when the classical velocity of the electron, v , approaches c, where relativistic effects can be dominant.In practical units, a 0 ¼ 0:86λ½μm ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I T ½10 18 W cm À2 p .In our case, a 0 = 0.16, which gives a relativistic correction of 0.6% to the electron mass in the THz field.Similar to the ponderomotive energy, the intensity required to reach the relativistic regime can be substantially reduced thanks to the Iλ 2 dependence.
Another important parameter relevant to ionization is the Keldysh parameter γ K ¼ ffiffiffiffiffiffi , where U i is the ionization potential energy of the atom or ion 45 .γ K can also be expressed as the square root of the ratio of the ionization time t ion to the THz period T. For a hydrogen atom (U i = 13.6 eV) at I T = 8.9 × 10 13 W cm −2 , γ K = 0.046.Such a small γ K value implies that ionization occurs so rapidly, compared to the THz period of 67 fs, that the THz field suppressing the Coulomb potential barrier is seen as a static field by the bound electron.Under this circumstance, the bound electron can tunnel through the suppressed barrier and be freed.In this regime of γ K < 1, tunneling ionization is a dominant mechanism over multiphoton ionization.

THz-induced ionization and damage
For the study of intense THz-matter interactions, we first generated THz waves by using a 25-μm-thick LN crystal doped with 5% MgO.For simplicity, the target sample was placed outside the vacuum chamber, as shown in Fig. 2. The energy of the THz waves emitted from the LN was estimated to be 0.62 mJ (with 2.3 J pumping) but drops to 0.27 mJ on the target after the waves pass through the ITO decoupler, the KBr window, and the Brewster-angled Si wafer.In addition, the pulse is chirped by the KBr window, providing GDD = − 1.15 × 10 5 fs 2 at 15 THz.This stretches the THz pulse duration from 0.23 ps to 1.3 ps.With E = 0.27 mJ, w = 37 μm, and τ = 1.3 ps, the peak intensity of the focused THz waves is estimated to be I T = 0.91 × 10 13 W cm −2 at the focus, and its corresponding THz field strength is E T = 83 MV cm −1 .The field strength is still high enough to induce tunneling ionization with the Keldysh parameter γ K = 0.2 less than 1.
As the first demonstration, a blank sample holder made of aluminum (Al) was exposed to a single THz pulse focused by the OAP mirror.Figure 7a shows a photograph of the aluminum target and a flash of plasma fluorescence captured by a charge-coupled device (CCD) camera (see Fig. S5 in Supplementary Material for raw photographs).Figure 7b shows a black aluminum foil exposed to an intense THz pulse and its resulting plasma emission.Figure 7c shows an InSb wafer together with its plasma emission.In these samples, the focused THz field is high enough to ionize the target by tunneling ionization and produce plasma on the target surface.For Al (U i = 6.0 eV) and In (U i = 5.8 eV), γ K < 0.1, confirming that tunneling ionization is dominant over multiphoton ionization.The electrons initially created by tunneling ionization are then heated via inverse bremsstrahlung absorption of the trailing THz waves over a picosecond time scale.Simultaneously, it undergoes avalanche ionization under strong THz fields, ultimately emitting plasma fluorescence by recombination, electronic transition, and bremsstrahlung processes beyond the picosecond time scale.Additionally, plasma expansion, ablation, and structural deformation occur on nanosecond time scales.
As the fourth target, a free-standing polymer (PET) film was exposed to intense THz radiation at the focus.This resulted in irreversible damage, such as the removal or destruction of the film, as shown in Fig. 7d.It shows clear holes drilled through the PET film after the exposure.Figure 7e, f show THz-induced damage created on the surfaces of GaAs and InSb wafers after 5 consecutive shots of THz radiation.To further investigate the dependence of THz-induced damage on the intensity of irradiated THz waves, the GaAs sample was exposed to 5 shots of THz pulses with the sample position varied along the THz beam direction (z) with an increment of Δz = 0.2 mm. Figure 8 shows a series of THz beam profiles (left) captured by the microbolometer and their corresponding images of THz-induced ablation (right) taken under a microscope.The damage patterns (right) closely resemble the intensity distributions of the irradiated THz radiation (left).The damage pattern, however, becomes less noticeable when the sample is positioned far away from the focal point, which drops the THz intensity and thus yields significantly less ionization.

Discussion
In conclusion, we demonstrated the generation and characterization of mJ-level, intense THz fields, approaching 260 MV cm −1 , at 15 THz from large-area phase-matched OR in LN by using a 150-TW Ti:sapphire laser.By testing various LN wafers with thicknesses of 25 μm, 50 μm, and 100 μm, also doped or undoped with 5% MgO, we found that the thinnest LN yields the largest THz energy owing to the linear absorption of THz and the cascade effect in LN.In addition, MgO-doped LNs produced 1.4 to 1.6 times more THz energy because of their higher nonlinearities than the undoped LNs.The spatial profiles of the produced THz waves, when focused by an OAP mirror, were characterized by a microbolometer camera, which measured the focal spot size as small as 43 μm in FWHM.The camera was also used in a single-shot THz interferometer, which allowed us to measure the field autocorrelation in time and its corresponding spectrum via a Fourier transformation.A numerical simulation was also conducted by solving 1-D coupled FMEs, and its results support the measured THz pulse spectrum centered at 14.9 THz with an FWHM of 1.9 THz.Furthermore, we demonstrated THz-driven tunneling ionization and subsequent collisional heating and avalanche ionization by focusing intense 15-THz waves onto aluminum, polymer and semiconductor samples.
Although the field strength of 260 MV cm −1 is the highest demonstrated thus far at 0.1 ~20 THz, there is still room for further enhancement.First, cryogenic cooling of the LN wafer can reduce the linear absorption of radiation at 15 THz, enabling the use of a thicker LN emitter to produce even higher THz energy.Second, our THz decoupler, based on ITO-coated windows and a Si Brewster-angled wafer, exhibits a throughput of 52%.A more efficient THz decoupler can deliver even higher THz energy onto the target.Third, tighter focusing can enhance the field strength even more.In our experiment, an OAP mirror with f = 4" was used to have space for

THz characterization
The energy of THz waves was measured by using a pyroelectric detector (Gentec, THz5D-MT-BNC).To block any residual laser beam leaking through the ITO decoupler, a 7.3-μm longpass filter was used in front of the pyroelectric detector.In addition, a 15.7-THz bandpass filter with a 2.4-THz FWHM was used to transmit THz waves within the band and cut off 0.1 ~3 THz waves also emitted from LN 36 .Such low-frequency THz waves were additionally cut off by the thick KBr window that strongly attenuates radiation at 40 − 300 μm 47 .One or more Si filters were placed in the beam path to attenuate THz energy and prevent any possible damage to the detector.
To characterize the profiles of the generated THz waves, we employed a microbolometer imaging sensor (camera).The camera was a commercial uncooled vanadium oxide (VOx) microbolometer (FLIR, Tau2 336), featuring a resolution of 336 × 256 pixels with the pixel pitch of 17 μm.It was connected to a Camera Link expansion board (FLIR) and a frame grabber (NI, PCIE-1433) to capture 14-bit images at 60 frames per second.This microbolometer camera was previously used to characterize THz waves over a wide range of THz frequencies 19,42,48,49 .
To acquire the temporal and spectral profiles of the THz waves, we implemented a single-shot THz spectrometer 42 .This consisted of a modified Mach-Zehnder interferometer and a microbolometer camera, as illustrated in Fig. 2. Two identical high-resistivity silicon wafers served as beam splitters (BS1 and BS2) to divide and combine the incident THz beams with equal dispersion and intensities.One of the beams passed through a fixed optical path, while the other passed through an adjustable optical path using a delay stage.Two identical OAP mirrors of 3-inch diameter and 9-inch focal length were used to separately focus the split THz beams onto the microbolometer sensor.To create an interference fringe pattern on the camera, the two beams were crossed at an angle on the sensor plane.The autocorrelation obtained from the interference pattern allowed us to characterize the THz pulse duration and obtain the spectral power of the waves via a Fourier transformation.
For further investigation of THz-induced ionization and damage, the generated THz waves were focused onto solid materials such as aluminum, polymer films, and semiconductor wafers by an OAP mirror with a focal length of 4 inches.To block any residual IR laser in the beam and simultaneously allow maximal THz transmission, a high-resistivity (>10 kΩ), thin (0.675-mm thickness) Si filter was placed in the beam path at a Brewster angle of 67°, allowing maximal THz transmission of 81% at 15 THz 34 .A similar Brewster-angled Si filter was used inside the vacuum chamber to completely block the laser after the ITO decoupler.

Simulation
We solved one-dimensional (1-D) coupled forward Maxwell equations (FMEs) for both THz and laser fields in the frequency domain 34,36 .This calculation fully resolves all carrier frequencies for both laser and terahertz pulses and is more powerful than nonlinear envelope equations commonly used to compute THz waves.Our 1-D coupled FMEs are given as where E represents the electric field, and the subscripts T and L denote THz and laser, respectively.Ω and ω are the angular frequencies of the THz and laser fields, respectively, and ζ ¼ z À v g t corresponds to the coordinate moving at the group velocity of the laser.The wavevector k is given by k ω ð Þ ¼ n ω ð Þω=c.On the right side of Eqs. ( 5) and ( 6), the first terms α T and α L and the second terms =c account for the absorption and dispersion of the THz and laser fields, respectively.The last term in Eq. ( 5) represents the OR process, a radiation source.The effective nonlinear coefficient d eff = d eff (Ω) is obtained from Ref. 35 .In Eqs. ( 5) and ( 6), THz absorption by free carriers, induced by threelaser-photon absorption in LN, is not considered because of its small effect (13 cm −1 ) compared to single-photon absorption at 15 THz in LN (~0.14 μm −1 ) 34,41 .The third and last terms in Eq. ( 6) represent the cascading effect and the self-phase modulation (SPM) of the laser, respectively.The cascading term represents a THz-induced modulation of the laser pulse.The SPM term involves the use of a nonlinear refractive index of n 2 = 10 −6 cm 2 GW −1 50 , with FT denoting a Fourier transformation.The input laser pulse was set to be Gaussian with a center wavelength of 800 nm and an FWHM of 78 nm.

Fig. 2
Fig. 2 Experimental setup for generating, characterizing, and focusing 15 THz waves to induce THz-driven ionization.The energy and spatial beam profile of the THz waves (blue beams) produced from a 3-inch-diameter LN wafer are measured with a pyroelectric detector and a microbolometer (middle dashed box).The temporal profile and spectral power are characterized by a single-shot THz interferometer (right dashed box).The generated THz waves are focused onto a solid target to demonstrate THz-induced ionization (left dashed box)

Fig. 3
Fig. 3 THz output energy scaling.a, b THz output energy measured as a function of input laser GDD at fixed laser energy of 2.3 J (except the green scatters taken at 3.4 J), obtained from x-cut LNs (a) with and (b) without 5% MgO doping.c, d THz output energy versus input laser energy obtained from x-cut LNs (c) with and (d) without 5% MgO doping at the optimal GDD for each LN.The insets show log-log plots.The gray lines represent a log-log slope of 2

Figure
Figure 5a-d present an example of THz beam profiles obtained in four distinct imaging modes.Figure 5a displays an intensity profile I x; y ð Þ produced when both THz beams were unblocked.It clearly shows an interference pattern with vertical fringes.Figure 5b, c show individual

Fig. 4
Fig. 4 THz beam profiling.a THz beam profiles captured at various positions along the propagation when focused by an OAP mirror.b THz beam profile at the focus and its lineouts and Gaussian fits in the horizontal (right top) and vertical (right bottom) directions and C 0 are the scaling coefficients that minimize the quasi-DC frequency in the Fourier transform of I bg ðx; yÞ.

Fig. 5
Fig. 5 THz pulse characterization in time and frequency.Microbolometer images obtained with (a) two THz replica beams turned on, (b) the first beam on only, (c) the second beam on only, and (d) both beams off.e Horizontal lineouts from each image, where I bg (x) represents the background signal,I bg ðxÞ¼ C 1 I 1 x ð Þ þ C 2 I 2 x ð Þ À C 0 I 0 x ð Þ. fTHz autocorrelation signals, I(τ) (black line) and background-free ΔI(τ) (red line), and (g) their corresponding THz spectra obtained by Fourier transformations Fig. 5 THz pulse characterization in time and frequency.Microbolometer images obtained with (a) two THz replica beams turned on, (b) the first beam on only, (c) the second beam on only, and (d) both beams off.e Horizontal lineouts from each image, where I bg (x) represents the background signal,I bg ðxÞ¼ C 1 I 1 x ð Þ þ C 2 I 2 x ð Þ À C 0 I 0 x ð Þ. fTHz autocorrelation signals, I(τ) (black line) and background-free ΔI(τ) (red line), and (g) their corresponding THz spectra obtained by Fourier transformations

Fig. 6
Fig. 6 Simulation of THz pulse generation in LN. a Refractive index n THz (top, black line), absorption coefficient α T (bottom, blue line), and effective nonlinear coefficient d eff (bottom, red line) of LN versus frequency.The dotted horizontal line (top) represents the group index of refraction n g at 800 nm.The arrows indicate phase-matched frequencies.b Calculation of THz spectral power (false color) with increasing laser propagation distance in LN. c Simulated THz energy (false color) plotted as a function of the initial GDD and the propagation distance in LN.The co-plot is the simulated laser-to-THz conversion efficiency with (red line) and without (yellow line) the cascading and SPM effects.The scatters are our experimental results with (cyan circles) and without (orange squares) MgO doping in LN. d Simulated THz waveform after 25-μm propagation in LN.The insets show the waveform (top) after bandpass filtering at 12-18 THz and its corresponding intensity envelope (bottom, red line).e Comparison between simulated (red line) and measured (black line) THz autocorrelation signals and (f) their Fourier spectra

Fig. 7 Fig. 8
Fig. 7 THz-induced ionization on solids.a-c Photographs of (a) aluminum holder, (b) black aluminum foil, and (c) InSb wafer targets irradiated by intense THz radiation, with the resulting plasma fluorescence captured by a CCD camera in each inset.d-f Microscope images of THz-induced damage made on (d) polymer (PET) film, (e) GaAs wafer, and (f) InSb wafer targets